! ******************************************************
! Grids and Particles Initialization
! ******************************************************

module Initialization
    use Input
    use Injector
    implicit none

contains

    subroutine setInput()
        implicit none

        call setControl()
        call setBoundaries()
        call setPhysics()
        call setOutput()
    end subroutine setInput

    subroutine setDimensionless()
        implicit none

        ! Characteristic Mass [g]
        M_unit = me
        n_unit = neCharacter
        if(neCharacter == 0.0d0) n_unit = M_unit*c**2/4/pi/e**2/L_unit**2
        v_Alfen = B_unit/sqrt(4*pi*neCharacter*mp)
        Omega_p = sqrt(4*pi*e**2*neCharacter/M_unit)
        Omega_g = B_unit*e/M_unit/c
        lambda_De = sqrt(kB*TeCharacter*11605/4/pi/neCharacter/e**2)
        SkinDepth = c/Omega_p !(inertial length)
        v_th = sqrt(kB*TeCharacter*11605/M_unit)
        gyroradius = v_th/Omega_g
        Cs = sqrt(Gamma*kB*TeCharacter*11605/M_unit)

        ! Characteristic Electric field [statV/cm] statV = g^1/2cm^1/2s^-1
        E_unit = 0.0d0

        if(B_unit == 0.0d0)then
            v_unit = c
            time_unit = 1/Omega_p
            B_unit = M_unit*c/time_unit/e
            E_unit = M_unit*v_unit/time_unit/e
        else
            UseConv = .false. ! can't use Poisson's Eq
            v_unit = v_Alfen
            time_unit = 1/Omega_g
            E_unit = M_unit*v_unit/time_unit/e
        end if
        if(L_unit == 0.0d0) L_unit = v_unit*time_unit
        p_unit = M_unit*c**2/4/pi/e**2/time_unit**2
        T_unit = M_unit*v_unit**2/kB/11605

        ! Simulation Parameters
        TotalTime = TotalTime/time_unit
        xMinBox = xMinBox/L_unit
        xMaxBox = xMaxBox/L_unit
        yMinBox = yMinBox/L_unit
        yMaxBox = yMaxBox/L_unit
        
        dt = 1.0d0/Nt_unit!sqrt(mp/me)/Nt_unit  ! Time Step
        dx = 1.0d0/Nx_unit  ! Cell Size in X
        dy = 1.0d0/Ny_unit  ! Cell Size in Y

        DumpPeriod = DumpPeriod/time_unit
        ! if(myid==7) write(*,*) xInjectMin,xInjectMax
    end subroutine setDimensionless

    subroutine initializeParallel()
        implicit none
        integer(kind=4) i,j

        RunTime = 0.0d0  ! Run Time Initialization
        DumpNumber = 0  ! Number of Diagnostic Datas

        if(UseConv .and. UseICS)then
            ncx = larger(INT(Nx_unit)*6,1+ceiling(lambda_De/L_unit*CxDebye*Nx_unit))
            ncy = larger(INT(Ny_unit)*6,1+ceiling(lambda_De/L_unit*CyDebye*Ny_unit))
        elseif(UseConv)then
            ncx = 1+ceiling(lambda_De/L_unit*CxDebye*Nx_unit)
            ncy = 1+ceiling(lambda_De/L_unit*CyDebye*Ny_unit)
        elseif(UseICS)then
            ncx = Nx_unit*5
            ncy = Ny_unit*5
        else
            ncx = ceiling(2*(1+v_th/v_unit)*dt*Nx_unit)   ! for injector
            ncy = ceiling(2*(1+v_th/v_unit)*dt*Ny_unit)
            ! ncx = 1
            ! ncy = 1
        end if

        ppi = mod(myid,ppx)
        ppj = myid/ppx

        if(UseInjector)then
            inx = 2*ncx+3+10*Nx_unit ! increasing the interacting time
            iny = 2*ncy+3+10*Ny_unit

            if(xMinInjector)then
                xMinBox = xMinBox-inx*dx
            end if
            if(xMaxInjector)then
                xMaxBox = xMaxBox+inx*dx
            end if
            if(yMinInjector)then
                yMinBox = yMinBox-iny*dy
            end if
            if(yMaxInjector)then
                yMaxBox = yMaxBox+iny*dy
            end if
        end if
        NxBox = ceiling((xMaxBox-xMinBox)*Nx_unit)  ! Cell Number in X
        NyBox = ceiling((yMaxBox-yMinBox)*Ny_unit)  ! Cell Number in Y
        ! Grid Size in Each Processor
      
        Nx = ceiling(real(NxBox,8)/ppx)
        Ny = ceiling(real(NyBox,8)/ppy)
        
        xMaxBox = xMinBox+NxBox*dx
        yMaxBox = yMinBox+NyBox*dy

        xMin = xMinBox+ppi*Nx*dx
        yMin = yMinBox+ppj*Ny*dy

        if(ppi == ppx-1) Nx = NxBox-ppi*Nx
        if(ppj == ppy-1) Ny = NyBox-ppj*Ny

        xMax = xMin+Nx*dx
        yMax = yMin+Ny*dy

        if(UseInjector)then
            if(xMinInjector .and. ppi==0)then
                if(Nx<=inx)then
                    write(*,*) "ERROR: injector region too big !"
                    write(*,*) "At ",myid
                    stop
                end if
                xInjectMin = xMinBox
                xInjectMax = xMinBox+(inx)*dx
                yInjectMin = yMinBox
                yInjectMax = yMaxBox
            end if
            if(xMaxInjector .and. ppi==ppx-1)then
                if(Nx<=inx)then
                    write(*,*) "ERROR: injector region too big !"
                    write(*,*) "At ",myid
                    stop
                end if
                xInjectMin = xMaxBox-(inx)*dx
                xInjectMax = xMaxBox
                yInjectMin = yMinBox
                yInjectMax = yMaxBox
            end if
            if(yMinInjector .and. ppj==0)then
                if(Ny<=iny)then
                    write(*,*) "ERROR: injector region too big !"
                    write(*,*) "At ",myid
                    stop
                end if
                xInjectMin = xMinBox
                xInjectMax = xMaxBox
                yInjectMin = yMinBox
                yInjectMax = yMinBox+(iny)*dy
            end if
            if(yMaxInjector .and. ppj==ppy-1)then
                if(Ny<=iny)then
                    write(*,*) "ERROR: injector region too big !"
                    write(*,*) "At ",myid
                    stop
                end if
                xInjectMin = xMinBox
                xInjectMax = xMaxBox
                yInjectMin = yMaxBox-(iny)*dy
                yInjectMax = yMaxBox
            end if
        end if
        
        ! Allocate Grid Variables

        if (ncx<(Nx-1) .and. ncy<(Ny-1)) then

            allocate(x(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(y(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Ex(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Ey(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Ez(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bx0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(By0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bz0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bx(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(By(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bz(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Exg(-ncx:Nx+ncx,-ncy:Ny+ncy))
            allocate(Eyg(-ncx:Nx+ncx,-ncy:Ny+ncy))
            allocate(Bxg(-ncx:Nx+ncx,-ncy:Ny+ncy))
            allocate(Byg(-ncx:Nx+ncx,-ncy:Ny+ncy))
            allocate(Bzg(-ncx:Nx+ncx,-ncy:Ny+ncy))
            allocate(Bwx(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bwy(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bwz(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bfx(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bfy(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(Bfz(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1))
            allocate(ux(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Fluid Velocity in X
            allocate(uy(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Fluid Velocity in Y
            allocate(uz(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Fluid Velocity in Z
            allocate(ne(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Number Density
            allocate(ni(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! ion Number Density
            allocate(ne_prev(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Number Density
            allocate(ne_temp(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Number Density
            allocate(Te(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Temperature
            allocate(pe(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Pressure
            allocate(ne0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Number Density
            allocate(Te0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Temperature
            allocate(pe0(-ncx-1:Nx+ncx+1,-ncy-1:Ny+ncy+1)) ! Electron Pressure

            do j=-ncy-1,Ny+ncy+1
                do i=-ncx-1,Nx+ncx+1
                    ! Mesh Coordinates
                    x(i,j) = xMin+i*dx
                    y(i,j) = yMin+j*dy
                end do
            end do
            
        else
            write(*,*)"Your entered convolution/guard space size exceed grid size !"
            write(*,*) ncx,ncy,Nx,Ny
            stop
        end if
        
        Ex = 0.0d0
        Ey = 0.0d0
        Ez = 0.0d0
        Exg = 0.0d0
        Eyg = 0.0d0
        Bx0 = 0.0d0
        By0 = 0.0d0
        Bz0 = 0.0d0
        Bx = 0.0d0
        By = 0.0d0
        Bz = 0.0d0
        Bxg = 0.0d0
        Byg = 0.0d0
        Bzg = 0.0d0

        Bwx = 0.0d0
        Bwy = 0.0d0
        Bwz = 0.0d0
        Bfx = 0.0d0
        Bfy = 0.0d0
        Bfz = 0.0d0

        ux = 0.0d0
        uy = 0.0d0
        uz = 0.0d0
        ne = 0.0d0
        ni = 0.0d0
        ne_prev = 0.0d0
        ne_temp = 0.0d0
        Te = 0.0d0
        pe = 0.0d0
        ne0 = 0.0d0
        Te0 = 0.0d0
        pe0 = 0.0d0

        A_avg = 0.0d0
        Z_avg = 0.0d0

        ! Setting Random Seed
        call random_seed()

        ! Initialization of PseudoParticle Linked List
        call createHeadOf(IonListHead)
        call createHeadOf(InjectorListHead)
        LocalIons = 0

        call createHeadOf(MPIListHeadLeft)
        call createHeadOf(MPIListHeadRight)
        call createHeadOf(MPIListHeadUp)
        call createHeadOf(MPIListHeadDown)
        call createHeadOf(MPIInjectorListHeadLeft)
        call createHeadOf(MPIInjectorListHeadRight)
        call createHeadOf(MPIInjectorListHeadUp)
        call createHeadOf(MPIInjectorListHeadDown)
        MPIIonsLeft = 0
        MPIIonsRight = 0
        MPIIonsUp = 0
        MPIIonsDown = 0

        ! If Periodic Condition is Mismatch
        if(xMinBoundary==2 .or. xMaxBoundary==2)then
            xMinBoundary = 2
            xMaxBoundary = 2
        end if
        if(yMinBoundary==2 .or. yMaxBoundary==2)then
            yMinBoundary = 2
            yMaxBoundary = 2
        end if
    end subroutine initializeParallel

    subroutine initializeParticles()
        implicit none
        type(PseudoParticle), pointer :: PrevIon, PrevInjector
        integer(kind=4) :: i, j, n
        logical :: CellIsEmpty, InjectorIsEmpty
    
        PrevIon => IonListHead
        PrevInjector => InjectorListHead
    
        do j=0,Ny-1
            do i=0,Nx-1
                if(UseInjector)then
                    if(x(i,j)<xInjectMax .and. x(i,j)>=xInjectMin .and.&
                    y(i,j)<yInjectMax .and. y(i,j)>=yInjectMin)then
                        do n=1,NInjector
                            InjectorIsEmpty = isInjectorEmpty(i,j,n)
                            if(.not. InjectorIsEmpty)then
                                call createInjector(i,j,n,PrevInjector,InjectorListHead)
                            end if
                        end do
                    else
                        do n=1,NSpecies
                            CellIsEmpty = isCellEmpty(i,j,n)
                            if(.not. CellIsEmpty)then
                                call createIons(i,j,n,PrevIon,IonListHead)
                                ! call isInInjector2(PrevIon,InjectorIsInside) ! just copy inner particle
                                ! if(InjectorIsInside) call copyIon(InjectorListHead,PrevIon)
                            end if
                        end do  
                    end if
                else
                    do n=1,NSpecies
                        CellIsEmpty = isCellEmpty(i,j,n)
                        if(.not. CellIsEmpty)then
                            call createIons(i,j,n,PrevIon,IonListHead)
                            ! call isInInjector2(PrevIon,InjectorIsInside) ! just copy inner particle
                            ! if(InjectorIsInside) call copyIon(InjectorListHead,PrevIon)
                        end if
                    end do
                end if
            end do
        end do
        if(UseInjector) call copyListA2B(InjectorListHead,IonListHead)
    end subroutine initializeParticles

    function isCellEmpty(i,j,Species)
        implicit none
        integer(kind=4) :: i, j, Species
        logical :: isCellEmpty
        type(FluidParticle) :: Plasma
        real(kind=8) :: Sum

        Sum = 0.0d0
        call setPlasmas(x(i,j),y(i,j),Species,Plasma)
        Sum = Sum+Plasma%ni
        call setPlasmas(x(i+1,j),y(i+1,j),Species,Plasma)
        Sum = Sum+Plasma%ni
        call setPlasmas(x(i,j+1),y(i,j+1),Species,Plasma)
        Sum = Sum+Plasma%ni
        call setPlasmas(x(i+1,j+1),y(i+1,j+1),Species,Plasma)
        Sum = Sum+Plasma%ni

        if(Sum==0.0d0)then
            ! Cell is empty
            isCellEmpty = .true.
        else
            ! Cell is not empty
            isCellEmpty = .false.
        endif
    end function isCellEmpty

    subroutine createIons(i,j,Species,PrevIon,ListHead)
        implicit none
        integer(kind=4) :: i, j, Species
        type(PseudoParticle), pointer :: PrevIon, ListHead
        type(PseudoParticle), pointer :: Ion
        type(FluidParticle) :: Plasma
        real(kind=8) :: ThermalVelocity(3)
        real(kind=8) :: Position(2)
        real(kind=8) :: xCenter, yCenter
        integer(kind=4) :: m

        xCenter = 0.5d0*(x(i,j)+x(i+1,j))
        yCenter = 0.5d0*(y(i,j)+y(i,j+1))
        call setPlasmas(xCenter,yCenter,Species,Plasma)

        do m=1,Plasma%ppc
            Position = getIonsPosition(i,j) ! Setup.F90
            call setPlasmas(Position(1),Position(2),Species,Plasma)

            if(Plasma%ni>0.0d0)then
                allocate(Ion)
                Ion%Z = Plasma%Z
                Ion%A = Plasma%A
                Ion%x = Position(1)
                Ion%y = Position(2)

                ThermalVelocity = getMaxwell(Plasma%A,Plasma%Ti) ! Setup.F90
                Ion%vx = Plasma%ux/v_unit+ThermalVelocity(1)*c/v_unit
                Ion%vy = Plasma%uy/v_unit+ThermalVelocity(2)*c/v_unit
                Ion%vz = Plasma%uz/v_unit+ThermalVelocity(3)*c/v_unit

                Ion%weight = Plasma%ni/Plasma%ppc/n_unit

                nullify(Ion%next)

                PrevIon%next => Ion
                PrevIon => PrevIon%next

                if(associated(ListHead,IonListHead)) LocalIons = LocalIons+1
            end if
        end do
    end subroutine createIons

    subroutine initializeGrids()
        implicit none
        integer(kind=4) i, j, n
        real(kind=8) :: Temp(3), neCell
        type(FluidParticle) :: Plasma
        real(kind=8) :: tempsum

        tempsum = 0.0d0

        do j=-1-ncy,Ny+ncx+1
            do i=-1-ncx,Nx+ncx+1
                if(UseField)then
                    ! Electric Field
                    call setEfield(x(i,j),y(i,j),Temp)  ! Input.F90
                    Ex(i,j) = Ex(i,j)+Temp(1)
                    Ey(i,j) = Ey(i,j)+Temp(2)
                    Ez(i,j) = Ez(i,j)+Temp(3)
  
                    ! Magnetic Field
                    call setBfield(x(i,j),y(i,j),Temp)  ! Input.F90
                    Bfx(i,j) = Temp(1)
                    Bfy(i,j) = Temp(2)
                    Bfz(i,j) = Temp(3)
                end if
            end do
        end do

        Bx = Bwx+Bfx
        By = Bwy+Bfy
        Bz = Bwz+Bfz
        Bxg = mapBfield2Grid(Bx)
        Byg = mapBfield2Grid(By)
        Bzg = mapBfield2Grid(Bz)
        Bx0 = Bx
        By0 = By
        Bz0 = Bz
        Exg = mapEfield2Grid(Ex)
        Eyg = mapEfield2Grid(Ey)

        do n=1,NSpecies
            do j=-1-ncy,Ny+ncx+1
                do i=-1-ncx,Nx+ncx+1
                    call setPlasmas(x(i,j),y(i,j),n,Plasma)
                    ux(i,j) = Plasma%ux/v_unit
                    uy(i,j) = Plasma%uy/v_unit
                    uz(i,j) = Plasma%uz/v_unit

                    neCell = Plasma%Z*Plasma%ni/n_unit
                    ne0(i,j) = ne0(i,j)+neCell
                    pe0(i,j) = pe0(i,j)+neCell*Plasma%Te/T_unit
                end do
            end do
            tempsum = tempsum + Plasma%ppc
            A_avg = A_avg + Plasma%ppc*Plasma%A
            Z_avg = Z_avg + Plasma%ppc*Plasma%Z
        end do
        A_avg = A_avg/tempsum
        Z_avg = Z_avg/tempsum

        ne = ne0
        ni = ne0/Z_avg
        ne_prev = ne
        ne_temp = ne
        pe = pe0
        Te0 = pe0/ne0
        Te = Te0
    end subroutine initializeGrids
    
    function larger(a,b)
        implicit none
        integer(kind=4) :: a, b, larger
        if(a>b)then
            larger = a
        else
            larger = b
        end if
    end function larger

end module Initialization
